############################################################
#Parametry skriptu

# (2) Nazev metody, (nepovinny parametr - vychozi hodnota je nazev skriptu)
P.name <- "Metoda Efektivnich Srazek"

# (4) Vstupni veliciny (povinny parametr - musi byt uvedeny veliciny, jednotky a casove intervaly)
P.input.variables <- c("precipitation", "temperature")
P.input.units     <- c("millimeter", "Celsius")
P.input.interval  <- c("day", "day")
P.input.columns   <- c("Srazky", "Teploty")
P.input.files     <- c("caslav.dat")

# (5) Vystupni veliciny (povinny parametr - musi byt uvedeny veliciny, jednotky a casove intervaly)
P.output.variables <- c("EP", "EP_klouzave")
P.output.units     <- c("", "")
P.output.interval  <- c("day", "day")
P.output.columns   <- c("EP", "EP2")
P.output.files     <- c("Vyhodnoceni_sucha.txt", "Vyhodnoceni_sucha.txt")

soubor = "caslav.dat"

zac_kon_sucha = function(zac_zaznamu, sucho) { #sucho: TRUE, FALSE, NA - z toho se zjišťují začátky a konce
  sucho_zac = NULL
  sucho_kon = NULL
  such = 0
  for (d in zac_zaznamu:dni) {
    if (!is.na(sucho[d])) {
      if (sucho[d] && ((!is.na(sucho[d - 1]) && (!sucho[d - 1])) || is.na(sucho[d - 1]))) {
          such = such + 1
          sucho_zac = c(sucho_zac, d)
      }
    }
    if (is.na(sucho[d]) && !is.na(sucho[d - 1]) && sucho[d - 1])
      sucho_kon = c(sucho_kon, d - 1)
    else if (!is.na(sucho[d - 1]) && !is.na(sucho[d]) && sucho[d - 1] && !sucho[d])
      sucho_kon = c(sucho_kon, d - 1)
  }
  return(list(sucho_zac, sucho_kon, such))
}



data = read.table(soubor, header=TRUE)

usek = 1:nrow(data)
sra = data$Srazky[usek]
tep = data$Teploty[usek]
dni = length(sra)


tep_kladna = vector("numeric", dni)
kalendar = as.Date(as.character(data$Datum[usek]), format = "%Y-%m-%d")


for (d in 1:dni)
  tep_kladna[d] = ifelse(tep[d] < 0, 0, tep[d])

pred_dni = 365
vaha = vector("numeric", pred_dni)
EP = DEP = SEP = PRN = EDI = vector("numeric", dni)

for (dv in 1:pred_dni) {
  vaha[dv] = 0
  for (d in 1:(pred_dni + 1 - dv))
    vaha[dv] = vaha[dv] + 1 / (pred_dni + 1 - d)
}

EP[1:pred_dni] = NA 
DEP[1:pred_dni] =  NA 
for (d in (pred_dni + 1):dni) {
  EP[d] = 0
  for (pd in 1:pred_dni)
    EP[d] = EP[d] + vaha[pd] * sra[d - pd]
}

kalendar_den = as.character(kalendar, format = "%m%d")
MEP = tapply(EP, kalendar_den, mean, na.rm = TRUE)

#pro MEP na začátku roku klouzavý průměr z hodnot na konci roku
klpr_dni = 5 #lichý
tmp_MEP = vector("numeric", length(MEP) + klpr_dni - 1)
tmp_MEP[(klpr_dni):length(tmp_MEP)] = MEP
tmp_MEP[1:(klpr_dni - 1)] = MEP[(length(MEP) - klpr_dni + 2):length(MEP)]
tmp_MEP = filter(tmp_MEP, rep(1/klpr_dni, klpr_dni), sides = 1)
okraj = (klpr_dni - 1) / 2
MEP_kl_tmp = tmp_MEP[(klpr_dni + okraj):length(tmp_MEP)]
MEP_kl_tmp = c(MEP_kl_tmp, tmp_MEP[klpr_dni:(klpr_dni + okraj - 1)])
MEP_kl = MEP
MEP_kl[1:length(MEP_kl)] = MEP_kl_tmp[1:length(MEP_kl_tmp)]

#pro EP koncové hodnoty jako NA
EP_kl = filter(EP, rep(1/klpr_dni, klpr_dni), sides = 2)
sd_EP = tapply(EP, kalendar_den, sd, na.rm = TRUE)

suma_pro_PRN = 0
for (pd in 1:pred_dni) {
  suma_pro_PRN = suma_pro_PRN + 1 / pd
}
for (d in (pred_dni + 1):dni) {
  DEP[d] = EP_kl[d] - MEP_kl[[as.character(kalendar[d], format = "%m%d")]]
  PRN[d] = DEP[d] / suma_pro_PRN
}
sd_PRN = tapply(PRN, kalendar_den, sd, na.rm = TRUE)
sd_DEP = tapply(DEP, kalendar_den, sd, na.rm = TRUE)


sucho_dle_SEP = vector("logical", dni)
sucho_dle_SEP =  NA
sucho_dle_EDI = vector("logical", dni)
sucho_dle_EDI =  NA 

for (d in (pred_dni + 1):dni) {
  SEP[d] = DEP[d] / sd_EP[[as.character(kalendar[d], format = "%m%d")]]
  sucho_dle_SEP[d] = ifelse(SEP[d] < 0, TRUE, FALSE)

  EDI[d] = PRN[d] / sd_PRN[[as.character(kalendar[d], format = "%m%d")]]
  sucho_dle_EDI[d] = ifelse(EDI[d] < -1, TRUE, FALSE)
}

sra_normal = tapply(sra, kalendar_den, mean, na.rm = TRUE)
sra_sd = tapply(sra, kalendar_den, sd, na.rm = TRUE)

tep_normal = tapply(tep_kladna, kalendar_den, mean, na.rm = TRUE)
tep_sd = tapply(tep_kladna, kalendar_den, sd, na.rm = TRUE)

c_pro_DI1 = 1
c_pro_DI2 = 0.1

typ_sucha = c("SEP", "EDI")
typu_such = length(typ_sucha)

# ind = list(names = typ_sucha)
such = vector("integer", typu_such)

for (ts in 1:typu_such) {
  tmp = zac_kon_sucha(pred_dni + 1, get(paste("sucho_dle_", typ_sucha[ts], sep="")))
  sucho_zac = tmp[[1]]
  sucho_kon = tmp[[2]]
  such[ts] = tmp[[3]]
  
  df = data.frame(sucho_zac)
  df = cbind(df, sucho_kon)
  
  CNS = ANES = vector("integer", such[ts])
  APD = PDst = Tst = DI1 = DI2 = PRN_sum = PRN_mean = PDst_mean = DI1_mean = DI2_mean = EDI_sum = EDI_mean = vector("numeric", such[ts])
  
  for (su in 1:such[ts]) {
    CNS[su] = sucho_kon[su] - sucho_zac[su] + 1
    ANES[su] = sum(SEP[sucho_zac[su]:sucho_kon[su]])    

    APD[su] = PDst[su] = Tst[su] = PRN_sum[su] = EDI_sum[su] = 0
    for (d in sucho_zac[su]:sucho_kon[su]) {
      kal_den = as.character(kalendar[d], format = "%m%d")
      APD[su] = APD[su] + sra[d] - sra_normal[[kal_den]]
      PDst[su] = PDst[su] + (sra[d] - sra_normal[[kal_den]]) /  sra_sd[[kal_den]]
      Tst[su] = Tst[su] + (tep[d] - tep_normal[[kal_den]]) /  tep_sd[[kal_den]]
      PRN_sum[su] = PRN_sum[su] + PRN[d]
      EDI_sum[su] = EDI_sum[su] + EDI[d]
    }
    DI1[su] = PDst[su] - c_pro_DI1 * Tst[su]
    DI2[su] = PDst[su] - c_pro_DI2 * Tst[su]
    PRN_mean[su] = PRN_sum[su] / CNS[su]
    PDst_mean[su] = PDst[su] / CNS[su]
    EDI_mean[su] = EDI_sum[su] / CNS[su]
    DI1_mean[su] = DI1[su] / CNS[su]
    DI2_mean[su] = DI2[su] / CNS[su]
    EDI_mean[su] = EDI_sum[su] / CNS[su]
  }
  df = cbind(df, CNS, ANES, APD, PDst, Tst, DI1, DI2, PRN_sum, PRN_mean, PDst_mean, DI1_mean, DI2_mean, EDI_sum, EDI_mean)
  
#   ind[[typ_sucha[ts]]] = df

  sloupcu = ncol(df)
  sink(paste("Charakteristika_obdobi_such_dle_",typ_sucha[ts], ".txt", sep=""))
  for (sl in 1:sloupcu) {
    cat(names(df)[sl])
    cat("\t")
  }
  cat("\n")
  for (su in 1:such[ts]) {
    for (sl in 1:sloupcu) {
      if (names(df)[sl] == "sucho_zac" || names(df)[sl] == "sucho_kon")
        cat(paste(kalendar[df[su, sl]], "\t", sep = ""))
      else
        cat(paste(df[su, sl], "\t", sep = ""))
    }
    cat("\n")
  }
  sink()
}

sink("MEP_rocni_prumery.txt")
cat(paste("Den v roce", "MEP", "MEP klouzavy","Prumerne srazky", "Smerodatna odch srazek", "Prumerna teplota", "Smerodatna odch teploty",  sep="\t")) 
cat("\n")
for (d in 1:length(MEP)) {
 cat(paste(names(MEP)[d],  MEP[d], MEP_kl[d],sra_normal[d], sra_sd[d], tep_normal[d], tep_sd[d], sep="\t"))
 cat("\n")
}
sink()

sink("Vyhodnoceni_sucha.txt")
cat(paste("Od", "EP", "EP2", "DEP", "SEP", "EDI", "PRN", "Sucho_dle_SEP", "Sucho_dle_EDI", sep="\t"))  
cat("\n")
for (d in 1:dni) {
  cat(paste(kalendar[d], EP[d], EP_kl[d], DEP[d], SEP[d], EDI[d], PRN[d], sucho_dle_SEP[d], sucho_dle_EDI[d], sep="\t"))
  cat("\n")
}
sink()







